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ABSTRACT 


A  progress  report  on  three  topics  that  were  studied 
during  the  third  quarter  of  Contract  AF19 (628) -5981  is  pre¬ 
sented.  Preliminary  results  on  spectral  estimates  of  seis¬ 
mic  events  are  presented  and  two  principal  conclusions  are 
suggested:  noise  prediction  is  probably  not  a  useful  means 
of  reducing  the  variability  of  spectral  estimates,  and,  for 
the  one  large  surface-focus  event  processed  to  date,  there 
is  a  significant  amount  of  signal  energy  above  4  Hz.  The 
problem  of  detecting  nuclear  explosions  in  the  presence  of 
large  natural  events  has  also  been  considered  during  this 
period.  This  problem  is  discussed  here,  with  the  emphasis  on 
the  possibility  of  "steering"  nulls  at  the  natural  events. 
Difficulties  resulting  from  the  sampled  nature  of  the  seismic 
records  and  from  coda  reverberations  are  also  discussed  in 
this  context.  Finally,  the  location  of  epicenters  by  beam¬ 
splitting  with  a  IAS A  is  discussed.  The  technique  is  de¬ 
scribed,  the  sources  of  error  are  analyzed,  and  data  from 
three  seismic  events  are  presented. 
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SECTION  I 


INTRODUCTION 


This  is  the  third  quarterly  rejport  on  Contract  AF19(628)- 
5981.  It  covers  the  period  from  September  15  to  December  15, 
1966.  Three  topics  are  discussed  in  this  report.  Section  II 
consists  of  presentation  of  preliminary  results  using  the 
spectral  estimation  techniques  discussed  in  the  last  quarterly 
report.  The  suppression  of  signals  from  large  earthquakes  for 
the  purpose  of  detecting  lower  magnitude  events  that  are  close 
in  both  time  and  space  to  the  large  event  is  discussed  in 
Section  III.  The  coherent  use  of  the  elements  of  LASA  to  lo¬ 
cate  the  epicenter  of  events  with  this  array  alone  is  consider¬ 
ed  in  Section  IV,  and  the  results  of  some  preliminary  calcu¬ 
lations  are  presented.  As  these  three  topics  are  rather  dis¬ 
joint,  introductory  and  concluding  remarks  concerning  each 
topic  appear  separately  in  the  appropriate  sections. 

All  three  discussions  are  status  reports  only,  for  much 
more  work  remains  to  be  done  on  each  of  the  topics  discussed. 

In  the  case  of  spectral  estimation,  only  one  event  has  been 
studied  to  date.  This  one  event,  however,  shows  substantial 
signal-to-noise  ratio  out  to  at  least  4  Hz.  The  study  of 
large  event  suppression  has  reached  the  point  of  discarding 
array  nulls  as  impractical,  and  is  now  progressing  to  a  study 
of  continental  sized  arrays.  Finally,  first  efforts  at  "beam¬ 
splitting"  for  precise  epicenter  determination  from  LASA  are 
in  their  infancy  but  certain  methods  of  refining  these  calcu¬ 
lations  are  clearly  indicated  and  will  be  pursued  in  the 
immediate  future. 
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SECTION  II 


SPECTRAL  ESTIMATES  OP  SEISMIC  EVENTS 


This  chapter  consists  of  a  progress  report  on  work 
relating  to  the  estimation  of  the  energy  density  spectra  of 
seismic  events.  The  theoretical  basis  for  this  work  was  pre¬ 
sented  in  Section  III  of  the  last  report[l] .  Poliowing  the 
preparation  of  that  report,  the  computer  programs  necessary 
to  implement  the  desired  calculations  were  developed.  To 
date  these  calculations  have  been  applied  to  one  seismic 
record  which  included  a  large  surface-focus  event  and  to 
three  additional  records  (from  the  same  seismometer)  consist¬ 
ing  only  of  noise.  Some  of  the  results  of  these  calculations 
are  presented  in  this  chapter. 

It  should  be  stressed  that  these  are  very  preliminary 
results,  the  calculations  being  based  on  a  single  seismic 
event  and  only  a  few  noise  records.  However,  these  calcula¬ 
tions  are  sufficient  to  tentatively  conclude  that  noise  pre¬ 
diction  will  not  be  a  useful  tool  in  spectral  estimation  of 
seismic  events  and  to  give  some  sample  results  on  seismic 
spectra . 

This  chapter  begins  with  a  brief  review  of  the  theoreti¬ 
cal  discussion  of  the  previous  report.  In  this  review  the 
equations  are  restated  in  the  finite-discrete  formulation  ap¬ 
propriate  for  the  digital  computer  implementation.  This  is 
followed  by  a  brief  discussion  of  some  of  the  programming  con¬ 
siderations  involved  in  carrying  out  these  calculations.  The 
next  section  presents  the  experimental  results  of  the  noise 
prediction  calculation.  This  is  followed  by  a  discussion  of 
the  results  of  the  spectral  calculations.  These  calculations 
involve  two  choices  of  time  windows  and  were  performed  on 
seismic  records  before  and  after  high-pass  filtering.  The 
final  section  consists  of  tentative  conclusions  and  recommen¬ 
dations  for  continued  work  on  this  project. 


2.1  FINITE-DISCRETE  FORMULATION  OF  APPROPRIATE  EQUATIONS 

2.1.1  Noise  Prediction 

In  the  previous  report  it  was  suggested  that  the  vari¬ 
ability  of  the  spectral  estimates  could  be  reduced  by  sub¬ 
tracting  from  the  observed  waveform  an  estimate  of  the  noise 
component  that  is  based  only  on  the  noise  preceding  the  seis¬ 
mic  event.  In  that  report  the  formulation  of  the  noise-pre¬ 
diction  operator  was  based  on  a  continuous  representation  of 
the  noise,  as  this  was  more  convenient  for  theoretical  pre¬ 
dictions.  In  this  section  the  finite-discrete  formulation  is 
briefly  stated.  The  theoretical  content  of  this  section  is 
well-known  in  several  contexts  (see,  for  example,  reference  2). 
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Let  the  samples  of  a  stationary  Gaussian  process  be  given 
by  n, ,  For  convenience  assume  that  the  rms  value  of  the  sam¬ 
ples  is  1.  Then,  writing  the  correlation  coefficients  as  Pk 


ninj  =  p|i-J 


(l) 


Considering  an  estimate  of  the  form 

M 


n  =  7  n  .h. 

a  -I  l,a 


a  =  0,1. 


,N-1 


(2) 


the  problem  is  to  find  the  h.  that  minimize  the  mean  square 
error,  given  by  efl2:  1’ 


e 


2 


a 


(n  -n  ) 
a  a 


(3) 


The  parameter  M  in  Equation  (2)  represents  how  much  of  the 
previous  noise  waveform  will  be  used  in  calculating  the  esti¬ 
mate.  In  principle,  the  quality  of  the  estimate  should  in¬ 
crease  with  M.  In  order  to  keep  computation  times  within  rea¬ 
sonable  limits,  some  compromise  must  be  made  in  choosing  M. 
Furthermore,  since  there  is  some  question  about  the  actual 
stationarity  of  the  noise,  a  value  of  M  that  is  too  large 
could,  in  fact,  detract  from  the  quality  of  the  estimate.  In 
the  calculations  discussed  below,  values  of  M  up  to  200  (which 
represents  10  seconds  at  the  sampling  rate  of  20  per  second 
used  in  these  calculations)  have  been  used.  The  parameter  N 
in  Equation  (2),  which  represents  the  number  of  samples  to  be 
predicted,  is  usually  chosen  to  equal  the  number  of  sample 
points  in  the  time  window  used  for  spectral  calculations. 
Typical  values  for  N  in  these  calculations  are  64  and  128. 
Expanding  Equation  (3),  in  terms  of  the  expression  for  the 
noise  estimate  given  by  Equation  (2),  and  using  Equation  (l) 
to  substitute  for  the  appropriate  expectations,  yields 


e 


2 

a 


(no  - 


^  n_ihi  " 


M 


j-i 


n-jhj,«) 


7-2?  „-?  ■? 


■  V  -  2J1  nan-i  Ji1P|i-j|hl(«hj.=. 


-  p0  -  2 


M  M 
j 


i  p..h.  +  y  y  pi..ih,  h. 

1  =  1  0+1  1)«  4=2  |i-j| 


O) 


8 


9h 


k,a 


-2p 


a+k 


+  2 


M 

l  P 


i=l 


k  =  1  ,M 


M 

P ot+k  -  p|k-i|hi,a 


k  =  1,M 


(5) 
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The  M  linear  equations  given  in  Equation  (5)  specify  the  M 
values  of  h.  a  that  will  yield  a  minimum  value  of  the  mean 
square  error i  Equation  (5)  may  be  conveniently  rewritten  in 
matrix  form: 


p0  pl* *  *  PM-1 

h, 

1,« 

P  a+1 

Pi  P o • • • 

h0 

2so 

— 

. 

. 

• 

• 

• 

• 

• 

• 

• 

p0 

hM,a_ 

Pa+M 

-  — 1 

— 

Assuming  the  covariance  matrix  is  not  singular,  the  solution 
for  the  h's  may  be  written  immediately  from  Equation  (6)  in 
terms  of  the  inverse  of  the  covariance  matrix. 


( — 1  •  •  • 

JG 

—  — 

-1 

_ 

P0  pl*  •  •  PM-1 

Pa  +  1 

i 

a 

s 

PM-1  p0 

p  a+M 

Substituting  Equation  (5)  into  Eauation  (4)  yields  a  rela¬ 
tively  simple  expression  for  the  minimum  mean  square  error. 


M  MM 

P°~  2Jl  P“+ihi>a  +  ill  XP|i-J  Ihi>ahj  ,a 


M 


P0  Pa+jhj,a 


M 


V 


P“+j+hj3“ 


(8) 


In  carrying  out  the  noise  prediction  with  actual  seis¬ 
mic  data,  a  long  sample  of  noise  (typically,  two  or  three 
minutes)  preceding  the  event  is  used  to  estimate  the  correla¬ 
tion  coefficients.  Rather  than  trying  to  invert  directly  the 
(M-dimensional)  covariance  matrix,  which  would  not  be  practi¬ 
cal  in  applications,  an  algorithm  due  to  Levinson  [3]  is  used. 
By  using  this  algorithm,  the  h ^  a  may  be  calculated  for  all  a 
of  interest.  Examples  of  the  resulting  values  of  ea  are  pre¬ 
sented  in  a  later  section. 


2.1.2  Spectral  Calculations 

Eollowing  the  notation  of  the  previous  report,  but  writ¬ 
ing  equations  in  discrete  form  rather  than  continuous,  yields 
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zi  '  ai[xi+ni] 


(9) 


where  the  xi  and  are  the  signal  and  noise  contributions  to 
the  observed  samples,  the  a^  form  the  time  window,  and  the  z ^ 
are  the  samples  to  be  transformed.  To  date,  only  rectangular 
time  windows  have  been  used,  and  scale  factors  have  been  chosen 
so  that  the  a^  may  be  written  as: 


ai  = 


l 


1 

/N 

0 


0  i  i 


N-l 


elsewhere 


where  N  is  the  length  of  the  time  window, 
based  on  the  z.  is  given  by 

-j2”kl 


(10) 


The  Fourier  series 


N 


Jk 


=  I 


i=l 


z.e 


N 


(ID 


and  the  spectral  estimates  are  simply  |  ZjJ  .  By  analogy  with 
the  previous  report  (or  by  direct  calculation)  the  expected 
contribution  of  the  noise  to  the  spectral  estimates,  Pk,  may 
be  written  as 

+(N-l)  |  ,  -j-2lTSk 

P,_  =  l  (1  -  J4L.)  p  e 

s=- (N-l) 


N 


k 


N 


N-l  .  . 

o  v  /  -1  s-.  tfirsk 

=  P°  sli  "  N  PsCOS' 


N 


(12) 


where,  as  before,  the  ps  are  the  correlation  coefficients  esti¬ 
mated  from  a  long  run  of  the  noise  preceding  the  event.  Equa¬ 
tion  (12)  is  simply  the  smoothed  periodogram  that  corresponds 
to  the  choice  of  time  window  given  in  Equation  (10).  Spectral 
calculations  will  be  presented  in  this  report  in  the  form  of 
plots  of  |Zk|2  and  Pk  vs.  k,  or  frequency. 


2.2  PROGRAMMING  CONSIDERATIONS 

Some  programming  was  necessary  in  order  to  carry  out  the 
calculations  reported  in  this  chapter.  Considerable  use  was 
made  of  subroutines  available  from  the  Teledyne  staff  in  writ¬ 
ing  these  programs.  The  subroutine  COOL,  which  calculates 
Fourier  series  by  means  of  the  Cooley-Tukey  algorithm  [4]  and 
the  subroutine  C00LV0LV,  which  convolves  two  waveforms  by 
calculating  and  multiplying  their  Fourier  series  and  then  cal¬ 
culating  the  corresponding  time  series,  were  used  extensively 
throughout  these  calculations.  (Descriptions  of  these  sub¬ 
routines  may  be  found  in  reference  .51 •)  For  the  noise  pre¬ 
diction  calculations  a  straightforward  implementation  of  the 
Levinson  algorithm  mentioned  above  was  used.  T’his  ■orograrr)  is 
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quite  similar  to  prediction  subroutines  available  from  Teledyne, 
but  it  is  sufficiently  different  that  it  was  easier  to  write 
it  separately  than  to  modify  an  existing  subroutine.  For  some 
of  the  calculations,  high-pass  filtering  of  the  seismic  record 
preceded  the  spectral  calculations.  To  carry  out  the  filtering, 
a  subroutine  was  written  which  calculates  the  (time-domain)  im¬ 
pulse  response  from  a  desired  frequency-domain  response.  Fil¬ 
tering  was  effected  by  using  COOLVOLY  to  convolve  the  impulse 
response  of  the  filter  with  the  entire  seismic  record.  For  the 
number  of  sample  points  used  in  these  calculations,  it  was  not 
necessary  to  segment  the  waveforms,  as  suggested  for  some  fil¬ 
tering  applications  by  Stockham  [6], 

One  problem  that  arose  in  filtering  the  noise  waveforms 
was  the  occurrence  of  spurious  zeroes  in  the  noise  waveform. 
When  parity  errors  occurred  in  the  digital  recording  of  the 
seismic  record,  a  value  of  zero  was  substituted.  These  spur¬ 
ious  zeroes  contributed  a  ringing  when  the  high  pass  filter 
was  used.  The  spurious  zeroes  were  removed  from  the  original 
record  by  the  following  simple  rule:  whenever  a  zero  or  a 
string  of  zeroes  occurred,  the  polarity  of  the  samples  imme¬ 
diately  preceding  and  following  the  zero(s)  was  checked;  if 
the  polarity  changed,  it  was  considered  a  "legitimate"  zero; 
if  the  polarity  did  not  change,  the  zero  was  replaced  by  the 
straight-line  interpolation  of  the  surrounding  values.  The 
noise-autocorrelation  function  was  estimated  by  using  C001V01V 
to  convolve  either  a  two-  or  a  three-minute  sample  of  the  noise 
with  itself.  Before  any  calculations  were  performed,  the  mean 
value  of  the  noise  was  removed  from  the  entire  record  and  the 
entire  record  was  divided  by  the  rms  level  of  the  noise.  Thus, 
for  all  of  the  calculations,  the  noise  process  had  a  zero  mean 
and  a  unity  mean  square  value. 


2.3  NOISE  PREDICTION  RESULTS 

The  noise  prediction  calculations  discussed  in  the  pre¬ 
vious  section  have  been  applied  to  the  seismic  record  contain¬ 
ing  the  surface  focus  event  and  also  to  three  other  seismic 
records  consisting  only  of  noise.  All  of  these  data  are  from 
the  same  seismometer,  which  is  the  center  element  of  LASA 
cluster  AO.  For  the  noise-prediction  calculations  reported 
here  the  autocorrelation  function  was  based  on  a  three-minute 
sample  of  the  noise  waveform  except  for  one  case  in  which  a 
two-minute  sample  was  used.  From  the  correlation  coefficients, 
the  hi, a  discussed  above  were  calculated.  The  hi  a  were  then 
used  to  calculate  the  predicted  noise  during  the  interval  of 
interest  and,  together  with  the  p^,  to  calculate  the  expected 
root-mean-square  root  of  Equation  (8)  above.  Figure  la  pre¬ 
sents  a  portion  of  the  autocorrelation  function  calculated 
from  the  noise  preceding  a  large  surface-focus  event.  The 
predicted  rms  error,  based  on  this  autocorrelation  function, 

Is  presented  in  Figure  2  along  with  the  corresponding  curve 
calculated  for  three  other  noise  samples.  The  bottom  curve  in 
Figure  2  is  the  one  corresponding  to  the  autocorrelation  func- 
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tion  of  Figure  la.  In  all  cases  the  noise  is  normalized  to 
have  an  rms  value  of  1,  so  that  these  curves  should  asymptote 
to  a  value  of  1.  For  the  bottom  curve  the  value  at  a  =  79, 
which  represents  a  "prediction  distance"  of  4  seconds,  is  0.85. 
At  one  second,  the  value  is  0.49;  at  two  seconds  it  is  0.73. 

The  other  three  curves  all  indicate  a  larger  rms  value  at  these 
times.  It  is  anticipated  that  the  time  windows  of  interest 
will  be  at  least  three  seconds  in  length.  From  Figure  2  it 
would  appear  that  the  noise  prediction  would  significantly  re¬ 
duce  the  rms  value  of  the  noise  during  only  a  small  fraction 
of  this  interval.  These  calculations  have  also  been  performed 
on  noise  waveforms  after  high-pass  filtering.  (The  details  of 
the  high-pass  filter  are  presented  in  the  next  section.)  The 
resulting  rms  error  curves  all  were  within  90$  of  their  asymp¬ 
totes  after  one  second. 

The  rms  error  curve  is  perhaps  the  simplest  measure  of 
the  quality  of  the  noise  prediction.  A  more  appropriate  mea¬ 
sure,  however,  is  the  reduction  in  the  variability  of  the 
spectral  estimate  that  can  be  achieved  by  using  the  noise 
prediction.  It  seems  unlikely  with  rms  errors  such  as  those 
of  Figure  2  that  any  significant  reduction  in  the  variability 
of  the  spectral  estimates  would  result  from  using  noise  pre¬ 
diction.  A  more  direct  measure  of  the  value  of  noise  pre¬ 
diction  in  this  context  would  be  a  theoretical  calculation, 
based  on  the  calculated  autocorrelation  function,  of  the  re¬ 
duction  in  the  variability  of  spectral  estimates  as  a  function 
of  frequency.  This  calculation  has  not  yet  been  performed 
with  these  data,  but  it  is  anticipated  that  this  calculation 
will  show  the  noise  prediction  not  to  be  of  much  value.  On 
physical  grounds  it  would  be  expected  that  the  noise  prediction 
would  reduce  the  variability  of  spectral  estimates  mostly  at 
the  low  frequencies.  This  expectation  is  confirmed  in  part 
by  the  rms  error  calculations  for  the  high-pass  filtered  noise. 
The  only  calculations  performed  to  data  that  directly  reflect 
on  this  question  are  calculations  of  the  spectrum  of  64  and  128 
point  noise  samples  with  and  without  noise  prediction.  As 
would  be  expected  from  theoretical  considerations,  the  results 
of  these  calculations  are  highly  variable.  A  quick  review  of 
these  calculations  suggests  that  the  noise  prediction  may  reduce 
the  variance  of  the  spectral  estimates  by  as  much  as  one-half 
at  low  frequencies  but  only  negligibly  at  higher  frequencies. 
Based  on  these  preliminary  calculations  it  is  anticipated  that 
noise  prediction  will  not  be  of  value  in  reducing  the  varia¬ 
bility  of  spectral  estimates. 


2.4  RESULTS  OF  SPECTRAL  CALCULATIONS 

The  seismic  event  used  in  these  calculations  is  a  large 
surface  focus  event  that  occurred  on  13  February  1966.  (PLE 
Card  No.  9;  Epicenter:  49. 8N,  78. IE;  Magnitude:  6.3;  Time: 
0  4  57.7  GMT).  The  recording  is  from  the  center  element  of 
1ASA  cluster  AO.  Before  sampling,  the  seismogram  signal  was 
passed  through  an  aliasing  filter  which  is  flat  out  to  5  Hz 
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and  provides  30  dB  of  attenuation  of  10  Hz.  A  plot  of  the 
seismic  waveform  is  presented  in  Figure  lb  and  the  time 
window  boundaries  are  indicated.  From  this  plot  it  is  clear 
that  this  is  a  very  high  signal -to-noise  ratio  event. 


Results  of  spectral  calculations  on  this  event  are  pre¬ 
sented  in  this  section  as  plots  of  I Zv> 2  vs.  k,  or  frequency. 
The  appropriate  smoothed  periodogram  (P^)  is  also  given  on 
the  same  plot  to  allow  comparisons  of  signal-to-noise  ratio 
at  different  frequencies.  In  viewing  these  data  it  should 
be  recalled  that  the  event  treated  here  is  a  very  high  signal- 
to-noise  ratio  event  so  that  general  conclusions  about  signal- 
to-noise  ratios  as  a  function  of  frequency  would  be  inappro¬ 
priate.  Two  time  windows  were  used  in  these  calculations: 

64  points  and  128  points.  (Integral  powers  of  two  were  chosen 
for  convenience  in  using  the  Fourier  series  algorithm,  but 
other  values  could  also  have  been  used.)  The  calculations 
based  on  the  64-point,  or  3.2  second,  time  window  are  pre¬ 
sented  in  Figure  3  and  those  based  on  the  128-point,  6.4 
second,  time  window  are  presented  in  Figure  4.  Because  of 
the  normalization  used  in  defining  the  time  window  (Equation 
10),  the  vertical  scales  on  these  two  figures  should  be  com¬ 
parable.  Specifically,  the  smoothed  periodogram  of  Figure  3 
should  be  somewhat  "smoother"  than  the  one  in  Figure  4»  and 
this  appears  to  be  the  case.  From  these  figures  it  is  clear 
that  both  signal  and  noise  spectra  fall  off  very  rapidly  with 
frequency.  A  principal  contributor  to  this  fall-off  is  the 
aliasing  filters  used  in  recording  these  data. 


These  figures  indicate  a  signal-to-noise  (power)  ratio 
on  the  order  of  100  at  frequencies  as  high  as  3.75  Hz.  Since 
there  is  very  much  more  energy  at  low  frequencies,  there  is 
the  possibility  that  this  apparent  energy  at  high  frequencies 
might  simply  be  a  result  of  the  sidelobes  in  the  Fourier  trans¬ 
form  of  the  rectangular  time  windows .  As  a  precaution  against 
this  possible  artifact,  the  seismic  records  were  high-pass 
filtered  and  the  calculations  repeated.  The  frequency  response, 
H(f),  of  the  high-pass  filter  is  given  by 
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Because  the  impulse  response  of  the  filter  is  symmetric  about 
the  origin  in  the  time  domain,  the  same  location  of  the  time 
window  was  used  with  and  without  filtering.  Spectral  calcu¬ 
lations  based  on  the  filtered  seismogram  are  presented  in 
Figures  5  and  6.  As  before,  the  smoothed  periodogram  was  ob¬ 
tained  by  first  estimating  the  autocorrelation  function  from 
a  long  run  of  (filtered)  noise  and  then  calculating  the 
appropriate  smoothed  periodogram  corresponding  to  the  time  win¬ 
dow.  In  Figures  5  and  6  there  is  apparently  some  signal  energy 
below  2.5  Hz,  despite  the  fact  that  the  frequency  response  of 
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FIGURE  3 

64-POINT  SPECTRAL  CALCULATION— UNFIITERED  RECORD 
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FIGURE  4 

128 -POINT  SPECTRAL  CALCULATION— FILTERED  RECORD 
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FIGURE  5 

64 -POINT  SPECTRAL  CALCULATION— FILTERED  RECORD 
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FIGURE  6 

128 -POINT  SPECTRAL  CALCULATION-- FILTERED  RECORD 
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the  filter  was  zero  in  this  region.  This  is  a  consequence 
of  the  sidelobe  effect  mentioned  above.  It  may  be  observed 
from  these  figures  that  the  good  signal-to-noise  ratio  at 
high  frequencies  is  apparently  still  present.  Comparison  of 
the  filtered  with  the  unfiltered  spectra  indicates  that  the 
spectral  levels  of  the  seismic  event  in  the  region  of  3*75 
to  5  Hz  are  not  much  affected  by  the  filtering  and  that,  there 
fore,  the  spectral  energy  at  these  frequencies  is  not  a  side- 
lobe  artifact  resulting  from  the  low  frequency  energy.  Based 
on  these  calculations  it  would  appear  that,  for  this  event, 
good  signal-to-noise  ratios  are  available  at  frequencies  as 
high  as  5  Hz.  Several  qualifications,  however,  should  be 
stated  with  this  tentative  conclusion.  All  of  the  calcula¬ 
tions  presented  here  are  based  on  a  single  seismic  event. 
Calculations  on  several  other  events  would,  of  course,  be 
necessary  before  any  inferences  about  typical  spectra  could 
be  made.  While  it  does  appear  that  the  high  frequency  energy 
is  not  a  sidelobe  artifact,  there  may  be  other  possible  arti- 
factual  contributors  that  have  not  yet  been  detected.  For 
example,  the  data  presented  in  these  figures  involve  a  very 
large  dynamic  range.  The  question  of  how  this  range  compares 
with  the  dynamic  range  available  in  the  recording  equipment 
has  not  been  thoroughly  explored. 

As  a  control  against  possible  errors  of  scale  factor, 
either  in  the  theoretical  formulations  or  in  the  programming, 
the  spectral  calculations  reported  above  for  the  seismic  event 
were  also  performed  on  a  control  interval  preceding  the  event, 
which  included  only  the  noise.  The  filtered  record  was  used 
for  the  examples  presented  in  Figures  7  and  8.  As  before, 
these  calculations  are  presented  alongside  the  corresponding 
smoothed  periodogram.  As  discussed  in  the  previous  report, 
the  smoothed  periodogram  represents  the  expected  value  of  the 
spectral  calculation,  and  the  variance  of  the  spectral  calcu¬ 
lation  at  a  given  point  is  at  least  as  large  as  the  square  of 
its  expected  value,  so  highly  variable  results  should  be  ex¬ 
pected  from  these  calculations.  Tfc.is  is  certainly  the  case  in 
Figures  7  and  8.  It  does  appear,  however,  that  the  smoothed 
periodogram  is  a  reasonable  curve  for  the  mean  value  of  the 
spectral  calculation  and,  therefore,  it  would  appear  that  no 
serious  scale  errors  were  involved  in  carrying  out  these  calcu 
lations . 

One  curious  similarity  between  Figures  6  and  8  is  the 
deep  notch  in  the  signal  spectrum  that  occurs  just  above 
3.75  Hz.  As  the  spectral  calculation  presented  in  Figure  8 
is  expected  to  be  highly  variable,  it  may  be  simply  coinci¬ 
dence  that  a  deep  minimum  in  that  calculation  occurs  at  the 
same  place  as  the  deep  minimum  in  Figure  6.  However,  there 
are  other  possible  explanations  for  this  occurrence  as  well. 
Reflections  from  a  subsurface  layer  near  the  seismometer 
could  lead  to  a  series  of  nulls  in  the  frequency  domain.  The 
seismogram  for  this  event,  which  was  presented  in  Figure  lb, 
suggests  secondary  arrivals,  for  example,  at  1  and  2  seconds 
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FIGURE  7 

64-POINT  SPECTRAL  CALCULATION  FOR  CONTROL  INTERVAL — FILTERED  RECORD 
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following  the  first  motion.  If  the  spectral  null  is  due  to 
local  conditions,  it  would  be  expected  that  it  should  also 
show  up  in  the  smoothed  periodogram  calculated  for  this  noise. 
(This  assumes  that  a  major  fraction  of  the  noise  energy  is 
distant  in  origin.)  A  slight  null  in  the  smoothed  periodo¬ 
gram  at  the  same  frequency  does  appear  in  Figure  8.  To  explore 
the  possibility  that  this  null  might  have  been  "smoothed  out," 
another  smoothed  periodogram  was  calculated  from  the  same  auto¬ 
correlation  function  with  four  times  the  resolution  of  the 
smoothed  periodogram  presented  in  this  figure.  The  results 
indicated  a  null  at  the  appropriate  frequency  at  which  the 
smoothed  periodogram  was  approximately  one-third  of  its  value 
at  the  surrounding  points .  This  result  may  not  be  statistically 
significant,  however,  since  the  spectral  calculation  involved 
a  "maximum  lag"  of  512  points  and  the  total  noise  record  used 
to  calculate  the  autocorrelation  function  was  only  2400  points. 
Similar  calculations  were  performed  on  other  noise  records 
from  the  same  seismometer  to  see  if  a  null  occurred  at  the 
same  frequency.  Although  deep  nulls  do  occur  in  some  cases  at 
nearby  frequencies,  none  of  them  occurred  at  exactly  this  fre¬ 
quency. 


2.5  CONCLUSIONS 

As  has  already  been  stressed,  the  results  in  this  report 
must  be  regarded  as  preliminary.  Further  calculations  on 
additional  seismic  events  will  be  necessary  before  any  specific 
conclusions  would  be  appropriate.  Now  that  the  major  program¬ 
ming  efforts  are  out  of  the  way,  however,  it  should  be  relative¬ 
ly  easy  to  perform  similar  calculations  on  several  other  seis¬ 
mic  events.  Tentative  conclusions  and  plans  for  continuing 
this  project  in  the  immediate  future  are  given  below. 


2.5.1  Noise  Prediction 

From  the  calculations  presented  in  this  report,  it  does 
not  appear  that  noise  prediction  will  be  a  useful  tool  in 
estimating  the  energy  density  spectrum  of  seismic  events. 

Eefore  this  possibility  is  discarded,  however,  a  few  additional 
calculations  will  be  performed.  If  these  calculations  (which 
require  little  further  programming)  yield  the  expected  results, 
no  further  work  on  noise  prediction  will  be  done  in  this  context. 


2.5.2  Parameter  Choices 

The  calculations  reported  above  all  involve  rectangular 
time  windows  and  only  two  lengths  of  these  time  windows.  The 
questions  of  what  is  the  best  shaped  time  window  and,  more 
importantly,  what  is  the  appropriate  duration,  have  yet  to  be 
resolved.  Additional  calculations  to  be  performed  on  other 
seismic  events  will  involve  the  so-called  harming  (raised- 
cosine)  window.  Eventually  the  question  of  optimal  window 
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lengths  will  have  to  be  studied.  The  choice  of  the  window 
length  will,  of  course,  have  to  be  a  compromise  between  the 
long  window  that  would  yield  good  resolution  in  the  frequency 
domain  and  a  short  window  that  would  presumably  increase  the 
si gnal-to -noise  ratio. 


2.5.3  Spectral  Smoothing  to  Improve  Signal-to-Noise  Ratio 

The  event  considered  in  this  report  is  atypical  in  that 
it  has  a  very  high  si gnal-to -noise  ratio.  In  general,  it 
would  be  expected  that  a  seismic  event  might  have  a  good  signal- 
to-noise  ratio  at  very  low  frequencies  but  not  at  all  frequen¬ 
cies.  When  this  is  the  case,  the  issue  will  arise  as  to  how 
best  to  smooth  the  spectral  calculations  over  frequency  in 
order  to  get  reliable  estimates  of  the  spectrum  of  the  event. 
This  question  was  mentioned  in  the  theoretical  discussion  of 
the  previous  report,  where  it  was  indicated  that  further  infor¬ 
mation  on  the  detailed  spectra  of  both  the  events  and  the  back¬ 
ground  noise  would  be  necessary  before  rules  for  smoothing 
could  be  developed.  It  is  hoped  that  some  progress  on  this 
topic  will  be  possible  as  the  results  of  additional  calcula¬ 
tions  become  available. 
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SECTION  III 


DETECTION  OE  NUCLEAR  EXPLOSIONS  IN  THE  PRESENCE  OP 
LARGE  NATURAL  EVENTS 


If  a  country  decided  to  violate  an  underground  test 
ban  treaty,  one  of  the  strategies  it  would  be  likely  to 
choose  is  selecting  an  area  of  high  seismicity  for  the  test. 
There  are  two  reasons  for  such  a  selection:  the  nuclear 
test  might  not  be  detected  at  all  if  its  signal  arrives  at 
the  monitoring  stations  simultaneously  with  a  large  earth¬ 
quake,  and,  if  it  is  detected  it  might  be  mistaken  for  an 
earthquake.  \ 

Large  events  produce  signals  of  considerable'' duration 
at  the  receiving  instruments.  A  magnitude  six  earthquake, 
for  example,  is  likely  to  last  many  minutes  before  its  power 
decreases  by  the  equivalent  of  one  magnitude  unit.  This 
slow  decay  rate  of  the  power  at  the  receiver  is  the  princi¬ 
pal  factor  that  limits  monitoring  capabilities .  Certain  array 
and  signal  processing  techniques  can  be  used  to  effectively 
reduce  the  interference  of  the  earthquake.  It  is  the  purpose 
of  this  chapter  to  describe  some  of  these  techniques,  and 
also  to  report  results  of  their  application  to  date. 

The  two  techniques  for  suppression  of  earthquake  signals 
that  are  considered  here  make  use  of  the  directional  charac¬ 
teristics  of  arrays.  In  one  case,  one  or  more  nulls  of  the 
array  pattern  are  steered  in  the  direction  of  signals  from 
the  interfering  earthquake.  The  limitations  that  these  nulls 
impose  on  searching  nearby  regions  for  test  explosions  will 
be  discussed  together  with  the  quality  of  nulls  which  can  be 
obtained.  In  the  second  case,  the  efficacy  of  steering  the 
main  lobe  to  search  for  test  explosions  in  the  vicinity  of  the 
earthquake  -  suppression  of  the  earthquake  signals  to  be 
accomplished  by  the  sidelobes  of  the  array  -  will  be  consider¬ 
ed.  In  both  of  these  approaches  four  considerations  are  of 
interest:  the  relative  magnitudes  of  the  two  events,  the 

relative  geometry  of  the  two,  the  relative  time  of  occurence 
of  the  two,  and  the  gometrical  distribution  of  receiving 
stations  about  the  epicenters  of  the  two  events. 

Before  proceeding,  it  should  be  noted  that  the  second 
approach  to  suppression  of  the  earthquake  signals  (which 
requires  a  narrow  beam  and  a  low  sidelobe  level  so  that  an 
earthquake  near  the  shot  epicenter  would  not  contribute  sig¬ 
nificantly  to  the  array  output  when  the  array  is  "steered"  to 
the  shot)  requires  very  large  arrays  for  any  chance  of  success. 
The  3  dB  beamwidth  of  LASA  Montana,  for  example,  is  on  the 
order  of  10°  at  teleseismic  ranges  -  far  too  large  for  the 
purpose  intended  here.  Thus,  continental-sized  arrays  would 
be  necessary  in  order  to  search  for  test  explosions  detonated 
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in  seismically  active  areas.  The  initial  problem  in  studying 
this  technique  is  one  of  obtaining  data  in  suitable  format. 

We  have  selected  events  and  stations,  and  requested  analog-to- 
digital  conversions  of  these  records  for  computer  processing; 
only  one  earthquake  has  been  processed  to  date. 

The  bulk  of  the  effort  during  this  period  was  devoted  to 
the  "null  solution"  to  the  problem,  the  remainder  of  this 
chapter  consists  of  a  discussion  of  this  approach.  We  shall, 
in  order,  consider  the  near  field  and  teleseismic  patterns  of 
dipoles  of  seismometers,  and  the  effect  of  a  combination  of 
such  dipoles  into  an  array  of  LASA  dimensions.  Also  discussed 
will  be  the  effect  of  time  sampling  rate  on  the  effectiveness 
of  a  nulling  array,  and  attempts  tc  remove  coda  reverberations 
which  represent  leakage  signals  frcm  energy  trapped  within  re¬ 
flecting  layers. 


3.1  NEAR  El ELD  RESPONSE  OE  SEISMIC  DIPOLES 

A  pair  of  seismometers  the  outputs  of  which  are  delayed 
according  to  some  rule  and  then  subtracted,  will  be  referred 
to  a  seismic  dipole.  Let  us  consider  such  a  dipole  as  is 
shown  in  Figure  9.  The  array  possesses  a  null  in  the  direc¬ 
tion  0=0  (broadside).  The  angle  ©  where  the  dipole  response 
is  maximum  is  given  by: 

9  =  sin"^(i  |)  X  =  fundamental  wave- 

^  d  length  (~14  km) 

We  can  see  that  a  seismic  dipole  of  length,  d,  of  approx¬ 
imately  60  km  would  possess  a  first  maximum  in  its  output  for 
an  angle,  9,  as  close  as  7  degrees  to  the  null.  This  high  slope 
of  the  pattern  can  be  of  great  help  in  the  problem  we  are  ana¬ 
lyzing. 

If  we  choose  a  delay  t  so  as  to  have  a  null  in  the  direc¬ 
tion  6  ,  the  angle  0  ,  (angle  of  maximum  output,  closest  to  0_\ 

will  b8  given  by  ml  0)> 

9  =  sin-1  (-  —  +  sin  9  ) 

m-^  2d  o' 

There  are  many  angles  where  the  conditions  for  maximum  output 
are  satisfied,  and  they  are  given  by  the  general  expression 

®m.  =  sin_1(-  Siti  £  +  sin  e0)  J=0,l,2.... 

t) 

The  angles  for  which  the  output  is  zero  (a  null  in  the  pattern) 
are  given  by 
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FIGURE  9 
SEISMIC  DIPOLE 


FIGURE  10 

SEISMIC  DIPOLE  OUTPUT 


0-1  >- 


en.  =  sln"1(:t  “a  +  5in  ec>  j=0,l,2.... 

t) 

The  angle  9  will  "be  a  function  of  range  and  azimuth. 

The  characteristics  of  the  pattern  are  not  greatly 
affected  by  the  introduction  of  a  delay,  the  null  is  just 
placed  at  the  corresponding  6  and  for  small  0O,  the  shape 
of  the  pattern  will  remain  ve£y  close  to  the  corresponding 
pattern  for  ©  =0  shown  in  Figure  1C. 


3.2  TELES El SMI C  RESPONSE  OP  A  SEISMIC  DIPOLE 

Consider  a  seismic  dipole  formed  with  two  of  LASA's 
elements.  We  shall  steer  a  null  to  a  certain  point  on  the 
surface  of  the  earth  and  compute  the  response  of  the  dipole 
for  events  coming  from  the  region  of  the  null.  Figure  11 
presents  a  schematic  diagram  of  the  configuration  under  con¬ 
sideration. 

In  order  to  steer  the  null  to  N,  we  have  to  introduce  a 
delay  tq  and  then  substract  the  two  signals. 


t  =  k.  ‘d‘Sin©^ 
o  Ao  o 

where  k.  =  i  [  ---e— — 1=  slope  of  the  travel- 

ao  V  I  degree J  ti£e  ourve 

d  [degrees]  =  length  of  the  dipole 
The  pattern  we  expect  is  shown  in  Figure  12. 

First  null: 

The  locus  of  the  first  null  is  obtained  by  setting  the 
relative  delay  between  arrivals  =  t  . 

.*.  d*k  •  sin0  =  d*k.^  •  sin©^ 

A  Ao  o 

(Equation  of  locus  5  in  Figure  12) 


The  first  maximum 

In  order  to  have  the  maximum  output,  the  difference 
between  the  actual  delay  in  arrival  and  the  "built  in"  delay 
tq  should  be  equal  to  4-T  where  T  is  the  period  of  the  domin¬ 
ant  frequency. 


d 


k.  sin© 
A 


sin© 

o 


T 

2 


In  a  similar  fashion  we  can  find  the  conditions  for  the  rest 
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FIGURE  11 

SEISMIC  DIPOLE  LOCATION  FOR  PATTERN  COMPUTATION 
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FIGURE  12 

EXPECTED  PATTERN  FOR  SEISMIC  DIPOLE 
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of  the  contours.  The  following  parameter  values  were  used 

T  =  1  sec  g. 
d  =  64  km  =  J~  degrees 

J.  XX  •  Xk 

*0  "  620 
AB„  =  320° 

The  resulting  dipole  response  is  plotted  in  Figure  13;  the 
corresponding  region  of  the  total  pattern  of  the  dipole  is 
indicated  in  Figure  14.  Observe  that  if  we  steer  a  null  at 
a  point  0  in  Figure  14  (broadside  to  the  dipole),  then  all 
points  lying  on  the  great  circle  that  joins  0  and  the  array 
center,  C  will  also  produce  a  zero  response.  On  the  other 
hand,  if  we  steer  the  null  to  a  point  M  (on  the  axis  of  the 
dipole),  the  great  circle  joining  M  and  C  will  show  alternat¬ 
ing  maxima  and  minima. 

The  radial  beamwidth  is  considerably  greater  than  the 
transverse  beamwidth  for  a  null  placed  at  M.  As  a  matter  of 
fact,  a  dipole  of  approximately  200  km  would  have  the  second 
null  at  about  20°  from  the  first  null. 

We  can  recognize  that  the  response  at  two  dipoles  such  as 
shown  in  Figure  15  is  going  to  show  regions  of  maximum  output 
such  as  indicated  by  the  shaded  areas  in  Figure  15. 

In  order  to  have  circular  symmetry  about  the  null,  N, 
d2  needs  to  be  larger  than  d-i.  A  distance  dp  on  the  order  of 
200  km  is  required. 

We  expect  that  having  the  freedom  to  choose  elements  from 
a  very  large  array,  we  can  obtain  a  pattern  that  has  a  central 
null  surrounded  by  a  region  of  high  response.  The  dipoles  are 
required  to  have  different  lengths  and  orientations .  Obviously 
the  width  of  the  region  of  high  response  depends  on  the  size 
of  the  dipoles. 


3.3  EFFECT  OF  TIME  SAMPLING  RATE  ON  EFFICIENCY  OF  NULLING 

In  order  to  be  able  to  steer  a  complete  null  (-  •  dB) , 
two  conditions  have  to  be  satisfied;  the  signals  i  the  ele¬ 
ments  have  to  be  identical  and  the  steering  delay  has  to  align 
both  signals  exactly.  The  first  condition  is  me^  fairly  well 
by  the  first  seconds  of  the  P-arrival  but  not  by  She  P-coda. 

It  was  hoped  however,  that  the  length  of  the  usa  i.e  interval 
might  be  increased  by  deconvolution  (or  some  oth  :  scheme  for 
removal  of  reverberation)  of  the  records. 


The  second  condition  cannot  be  met  using  discrete  time 
records  without  interpolation,  since  the  desired  time  resolu- 
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LATITUDE 


Shaded  regions  =  All  points  within  3  dB  of  maximum 
d  =  6*1  km 

Location:  LASA  center 
Null  coordinates:  a  =  62° 

AZ  =  320° 
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FIGURE  13 
DIPOLE  PATTERN 


FIGURE  1H 

REGION  OF  THE  TOTAL  PATTERN  OF  THE  DIPOLE 
ILLUSTRATED  IN  FIGURE  13 


PATTERN  OF  AN  ARRAY  OF  TWO  DIPOLES 
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tion  is  less  than  the  sampling  period.  Using  a  monochromatic 
wave  of  1  Hz  for  illustration,  it  can  be  shown  that  the  rela¬ 
tive  suppression  will  be  given  by: 

G  =Y2 IT  cos  2tta'  where  A  misalignment 
in  seconds 


Three  sample  values  of  suppression,  corresponding  to  three 
possible  values  of  a  are  tabulated  below 


A 


Sample  pointl 
Lfor  20  s.p.s. 


G 


[seismic  . 
^magnitude 


0 

0 

0 

.00 

.00 

* 

.025 

.155 

-16 

-.5 

1 

.05 

.285 

-10.5 

-.5 

From  this  table 

it  is  clear 

that 

the  obtainable 

suppression 

is  seriously  deteriorated  by  timing  errors.  In  order  to  im¬ 
prove  the  suppression  capabilities  a  higher  sampling  rate  or 
an  interpolation  between  the  sample  points  is  required. 


3.4  IMPROVEMENT  IN  CODA  CORRELATION 

As  mentioned  earlier  in  this  chapter,  there  is  a  very 
high  correlation  between  the  first  few  seconds  of  the  P-wave 
signals  received  at  different  receiving  sites,  even  though 
these  sites  may  be  separated  by  large  distances.  It  is  likely, 
then,  that  a  high  degree  of  cancellation  of  such  signals  could 
be  achieved  by  subtracting  one  from  another  after  normalization 
of  the  peak  values  of  these  signals.  P-coda,  on  the  other 
hand,  is  only  partly  correlated  between  receiving  sites  for 
signals  of  shallow  earthquake  origin  (p»=0.4)  [/]  .  For  shots, 
the  average  coda  correlation  is  close  to  zero,  provided  the 
receiver  sites  are  separated  by  at  least  100  km.  It  is  hypoth¬ 
esized  that  there  is  a  correlated  part  of  the  earthquake  coda 
resulting  from  such  things  as  closely  spaced  after-shocks,  etc. 

Since  subtraction  of  uncorrelated  signals  will  not  produce 
the  desired  result  of  lowering  the  coda  energy,  attention  was 
turned  towards  reducing  that  part  of  the  coda  generated  in  the 
vicinity  of  the  receiver;  the  hope  was  that  the  resulting  coda 
signal  correlation  would  increase  to  such  an  extent  that 
"nulling"  would  be  a  useful  tool.  As  a  first  step  in  examina¬ 
tion  of  this  problem,  results  of  Texas  Instruments  application 

If  "deconvolution"  to  several  teleseismic  events  were  studied 
8J  .  This  study  indicated  that  the  results  of  these  calcula- 
ions  would  not  be  satisfactory  for  our  purposes.  It  was  con¬ 
cluded  that  the  difficulties  are  a  consequence  of  the  fact  that 
the  two-way  travel  time  in  any  trapped  layer  is  substantially 
less  than  the  duration  of  the  elementary  wavelet  of  the  seismic 
signal. 
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FIGURE  16 

MODEL  FOR  RECEIVING  SITE 


FIGURE  17 

EQUIVALENT  FILTER  FOR  THE  LAYERED  EARTH 
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(POSITIVE  REFLECTION  COEFFICIENT) 


FIGURE  18 

RELATIONS  BETWEEN  DELAYS  IN  THE  TIME  FUNCTION 
AND  MINIMA  IN  THE  ENERGY  SPECTRA 
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In  summary,  it  can  be  stated  that  small  dipoles  (l  or  2 
wavelengths)  are  not  useful  for  the  purpose  of  separating  shots 
and  quakes,  since  the  nulls  produced  by  such  dipoles  cover 
practically  the  entire  surface  of  the  earth.  Further,  only  a 
small  number  of  acceptable  dipoles  can  be  formed  from  an  array 
such  as  1ASA  Montana.  Finally,  an  interpolation  of  seismic 
records  between  sample  times  would  be  necessary  if  nulling 
arrays  are  to  be  used  to  solve  the  overlapping  explosion- 
earthquake  problem. 

An  effort  was  made  to  surmount  this  problem  by  choosing  a 
different  calculation  for  the  removal  of  reverberation  in  the 
vicinity  of  the  receiver.  Here,  the  crust  was  assumed  to  be 
made  up  of  plane-parallel  layers,  vrhere  the  trapped  signals 
making  up  the  reverberation  are  as  shown  in  Figure  16.  It 
was  next  assumed  that  the  transmission  of  signals  through  this 
crustal  model  could  be  represented  by  a  cascade  of  feedback 
loops  as  shown  in  Figure  17-  The  delays  shown  in  this  figure 
are  equal  to  the  equivalent  of  the  two  way  travel  time  within 
layers,  and  the  loop  gains  are  taken  to  be  the  product  of  the 
reflection  coefficients  at  the  top  and  bottom  of  the  layers. 

If  the  assumed  propagation  model  were  correct,  the  receiv¬ 
ed  signal  y(t)  would  contain  sufficient  information  about  the 
t's  to  allow  identification  of  the  delays  and  gains  associated 
with  the  reflecting  layers.  Knowledge  of  these  parameters 
would  then  permit  construction  of  the  inverse  filter  for  the 
assumed  crustal  model.  This  inverse  filter  would  be  a  cascade 
of  feed- forward  loops  with  gains  and  loop  delays  identical  with 
those  representing  the  crust. 


3.4.1  Determination  of  the  layer  Delay  Times 

Where  a  time  function  is  made  up  of  repeated  versions  of 
some  wavelet,  with  even  spacings  in  time,  the  spectrum  of  the 
time  function  will  contain  nulls  whose  location  in  frequency 
is  related  to  the  time  spacing  of  the  wavelets  making  up  the 
time  function.  The  spacing  between  these  nulls  in  the  spectrum 
is  equal  to  the  reciprocal  of  the  time  separation  of  the  wave¬ 
lets.  The  frequency  of  the  lowest  null  is  dependent  on  whether 
the  time  function  consists  of  wavelets  with  all  the  same  polar¬ 
ity  (positive  reflection  coefficients)  or  consists  of  wavelets 
with  alternating  polarity  (negative  reflection  coefficients). 
These  qualitative  observations  are  illustrated  in  Figure  18. 

Of  course,  with  actual  events,  one  expects  several  reverbera¬ 
ting  layers.  These  would  produce  several  sets  of  nulls  in  the 
spectrum  of  the  received  waveform,  and  the  various  sets  would 
have  to  be  sorted.  The  sorting  problem  will  be  discussed  in 
more  detail  below. 
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3.4.2  Determination  of  Reflection  Coefficients 

Provided  it  is  possible  to  determine  the  delay  times 
associated  with  the  various  reverberating  layers,  the  re¬ 
flection  coefficients,  or  loop  gains,  may  then  be  determined 
by  passing  the  time  function  through  a  feed-forward  loop  with 
time  delay  equal  to  the  longest  observed  layer  delay.  The 
gain  of  this  loop  is  varied  between  zero  and  one,  and  the 
energy  in  the  output  observed  during  a  time  interval  that  be¬ 
gins  after  the  onset  of  the  P-wave.  The  time  between  the  on¬ 
set  of  the  P-wave  and  the  beginning  of  this  time  interval  is 
chosen  to  be  the  larger  of  the  following  two  quantities:  the 
duration  of  the  wavelet  and  the  delay  time  of  the  loop.  Proper 
loop  gain  is  that  which  minimizes  this  energy.  The  original 
time  function  is  passed  through  the  loop  with  this  gain,  and 
the  output  time  function  is  then  passed  through  a  feed-forward 
loop  with  the  second  longest  delay  time.  The  process  is  re¬ 
peated  until  all  of  the  delays  inferred  from  the  spectrum  have 
been  used.  The  final  time  function  should  then  have  local  re¬ 
verberations  removed,  and  the  correlation  of  coda  for  two  sta¬ 
tions  should  be  considerably  larger  than  at  the  start  of  the 
processing. 

This  processing  scheme  was  checked  with  a  synthetic  seis¬ 
mogram.  Recovery  of  the  pre-crust  waveform  was  satisfactorily 
achieved.  The  next  step  was  application  of  the  technique  to 
the  record  of  an  actual  seismic  event. 


3.4.3  Test  of  the  Technique  for  Removing  local  Reverberations 
with  an  Actual  Event 

To  test  the  computational  technique  described  above  on 
actual  data,  a  surface-focus  event  (13  Peb.  1966,  PDE  Card  No.  9; 
Epicenter:  49. 8N,  78. IE;  Magnitude:  6.3;  Time:  0  4  57.7  GMT) 
recorded  at  LASA  was  chosen.  This  choice  was  predicated  on  a 
number  of  factors.  First,  the  event  was  one  of  large  magnitude. 
Thus,  the  signal-to-noise  ratio  was  large.  Second,  as  a  sur¬ 
face  focus  event,  the  seismogram  was  relatively  simple,  and 
the  P-coda  should  be  generated  almost  entirely  in  the  vicinity 
of  the  receiver.  Finally,  the  seismogram  shows  an  obvious  re¬ 
flection  approximately  one  second  after  onset  of  the  P-wave. 

The  energy  density  spectrum  (based  on  a  60  sec  sample  be¬ 
ginning  at  the  first  motion)  of  this  event  did  show  many  minima 
and  their  interpretation  in  terms  of  a  layer  structure  was 
attempted.  The  first  task  was  to  find  the  delays  t1  that  are 
consistent  with  the  minima  found.  For  this  purpose  a  program 
was  written  which  computes  the  expected  minima  (for  positive 
and  negative  reflections)  for  t's  between  1  sec  and  50  sec. 

These  t's  were  formed  by  a  geometric  series  of  ratio  1.02.  The 
number  of  matches  with  the  experimental  minima  was  determined 
and  three  values  of  delays  which  showed  a  high  correspondence 
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were  selected:  5.6  sec,  1.24  sec,  1.06  sec.  These  values 
corresponded  to  assumed  positive  reflections.  All  attempts 
with  feed-forward  loops  to  determine  the  reflection  coefficients 
of  the  three  layers  met  with  failure,  as  no  minimum  energy  was 
obtained  for  any  gain  within  the  extremes  of  one  and  zero. 
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SECTION  IV 


LOCATION  OE  EPICENTERS  BY  BEAMSPLITTING  WITH  A  LASA 


Traditional  seismology  provides  a  technique  for  the 
approximate  determination  of  the  epicenter  of  a  large  event 
from  the  records  available  at  a  ^single  seismic  station.  This 
method,  first  used  by  Galitzin  [9]  ,  entails  determination  of 
the  station-epicenter  bearing  from  the  P-wave  onset  signals 
from  a  three  component  seismometer,  and  the  determination  of 
station-epicenter  range  from  the  differential  arrival  times  of 
the  P  and  S  phases  from  the  event.  Unfortunately,  this  tech¬ 
nique  was  of  limited  use,  since  the  S  phase  of  a  teleseism 
cannot  be  identified  for  events  much  smaller  than  magnitude 
six.  Moreover,  the  precision  of  bearing  determination  with  a 
three  component  seismometer  is  on  the  order  of  several  degrees. 
The  advent  of  large  aperture  seismic  arrays  offers  the  possibi¬ 
lity  of  overcoming  both  of  the  shortcomings  of  the  Galitzin 
method,  and  in  fact,  making  possible  relatively  precise  epi- 
centeral  determinations  from  a  single  seismic  station. 

The  way  in  which  a  LASA  makes  this  possible  is  through 
determination  of  the  direction  of  arrival  of  the  phase  front 
of  arriving  seismic  signals,  the  P-wave  in  particular.  With 
an  apperture  some  ten  wavelengths  in  dimension,  rather  precise 
determinations  of  both  the  station-  epicenter  bearing  and  the 
angle  of  incidence  (directly  related  to  station-epicenter 
range)  can  be  made.  Here  it  is  important  to  point  out  the  dif¬ 
ference  between  the  resolution  and  the  precision  of  angular 
measurement  possible  with  an  array.  The  resolution  of  an  array, 
the  angular  separation  between  two  sources  required  before  the 
array  can  identify  the  fact  of  multiple  sources,  is  limited  en¬ 
tirely  by  the  aperture  of  the  array.  It  is  generally  taken  to 
be  approximately  equal  to  the  beamwidth  of  the  array.  This 
beamwidth  is  approximately  10°  for  an  array  of  the  dimensions 
of  LASA  Montana.  The  precision  of  angular  determination  possi¬ 
ble  with  an  array,  on  the  other  hand,  is  related  to  both  the 
array  aperture  and  the  signal-to-noise  ratio  for  the  received 
signals.  In  theory,  an  infinite  S/N  with  the  smallest  of  arrays 
would  permit  infinitely  fine  angular  measurements. 

Clearly  the  degree  of  precision  achievable  with  a  given 
aperture  and  a  given  S/N  will  depend  on  the  illumination  across 
the  aperture;  however,  this  effect  is  weak.  For  filled  arrays 
(receiving  elements  at  least  every  half  wavelength  throughout 
the  aperture)  and  continuous  wave,  analog  signals,  it  can  be 
shown  that  where  noise  alone  is  the  interference,  the  mean 
error  in  angular  determination  is  zero,  and  the  standard  de¬ 
viation  of  this  measurement  is  inversely  proportional  to  the 
square  root  of  the  power  signal  to  noise  ratio  at  the  output  of 
the  array  [lo]  .  While  there  are  other  sources  of  interference 
in  the  case  of  a  LASA  making  epicentral  estimates,  let  us  first 
proceed  to  discuss  the  implementation  of  a  computational  scheme 
to  achieve  angular  measurements  with  precision  near  the  limit 
just  discussed. 
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4-1  THE  MONOPULSE  TECHNIQUE  OE  BEAMS  PUTTING 


Consider  the  output  of  an  array  as  the  array  is  aimed 
successively  at  various  angles  when  there  exists  a  signal 
arriving  from  one  direction.  This  output  would  peak  when  the 
array  is  pointed  directly  at  the  incoming  signal,  and  would 
fall  off  as  the  array  was  steered  away  from  the  signal.  This 
array  pattern  is  illustrated  in  Figure  19.  In  the  absence  of 
noise,  one  could  ideally  steer  the  array  in  successively  finer 
steps  in  the  vicinity  of  the  maximum  response,  and  eventually 
determine  that  angle  of  the  signal  with  any  degree  of  precision 
desired.  When  noise  is  added  to  the  output  of  the  array,  how¬ 
ever,  the  point  of  maximum  response  to  the  signal  is  increas¬ 
ingly  obscured.  Thus  other  means  must  be  sought  to  determine 
the  angle  of  arrival  of  the  signal. 

This  problem  is  one  which  has  been  successfully  attacked 
in  the  case  of  tracking  radar  systems,  and  the  technique  devel¬ 
oped  is  directly  applicable  to  the  seismic  problem.  In  radar, 
the  improvement  in  angular  accuracy  over  conventional  use  of 
arrays  has  been  called  beamsplitting,  and  the  particular  tech¬ 
nique  of  beamsplitting  to  be  discussed  here  is  called  monopulse. 
The  monopulse  technique  makes  use  of  two  beams  aimed  a  small 
angle  apart.  The  output  of  one  of  these  beams  is  then  sub¬ 
tracted  from  that  of  the  other.  If  the  two  beams  are  identical 
in  shape,  this  output  becomes  zero  as  the  pair  of  beams  is 
scanned  when  the  signal  is  arriving  at  an  angle  midway  between 
the  two  beams.  Here  the  slope  of  the  array  pattern,  rather 
than  being  zero  as  it  is  for  the  single  array  being  scanned 
past  the  signal,  has  a  relatively  steep  slope.  Thus  it  is  far 
easier  to  implement  automatic  detection  of  the  zero  crossing  of 
the  array  output,  and  hence  the  angle  of  arrival  of  the  signal. 


4.2  SOURCES  OE  INTEREERENCE  SPECIFIC  TO  THE  LASA  PROBLEM 

In  addition  to  the  noise  problem  discussed  above,  there 
are  several  additional  aspects  of  the  seismic  problem  which 
will  introduce  errors  in  the  determination  of  angles  by  the 
monopulse  technique.  Most  important  of  these  is  the  fact  that 
there  exist  local  travel-tim^  anomalies  within  the  aperture  of 
the  array.  At  LASA  these  anomalies  have  been  found  to  be 
functions  of  both  the  bearing  from  which  the  signals  arrive  as 
well  as  the  range  between  LASA  and  the  epicenters  of  the  events. 
These  anomalies  have  been  found  to  be  as  great  as  several  tenths 
of  a  second  [ll]  .  This  magnitude  of  potential  time  errors  makes 
it  quite  difficult  to  form  beams  at  all,  let  alone  to  take  pre¬ 
cise  measurements  of  arrival  angles  from  the  difference  of  two 
beams.  Two  separate  approaches  are  being  taken  to  attempt  re¬ 
moval  of  these  travel-time  anomalies  from  records  for  processing. 
First,  where  the  S/N  is  large  enough,  these  anomalies  are  mea¬ 
sured  from  the  seismic  records  themselves,  and  second,  where 
this  is  not  possible,  the  average  of  anomalies  from  events  pre¬ 
viously  measured  for  the  particular  area  are  used.  For  the 


35 


I 

ARRAY 

OUTPUT 


FIGURE  19 

ARRAY  RESPONSE  TO  A  SINGLE  INPUT  SIGNAL 


36 


moment,  both  methods  are  being  used  for  large  S/N  events  and 
the  results  compared.  Not  enough  statistical  evidence  is 
available  at  this  time  to  predict  the  magnitude  of  errors 
likely  to  result  from  these  travel-time  anomalies. 

In  addition  to  the  travel-time  anomalies  discussed  above, 
the  time-sampled  nature  of  1ASA  signals  introduces  a  form  of 
quantization  noise  in  the  beam  patterns  calculated  from  actual 
seismic  records.  While  no  new  information  will  be  gained  from 
the  process,  it  is  planned  to  interpolate  the  seismograms  so 
as  to  have  a  sample  rate  of  100  per  second  internal  to  the 
computations.  During  the  period  covered  by  this  report,  how¬ 
ever,  this  interpolation  has  not  yet  been  included  in  the 
computations . 

Another  problem  associated  with  1ASA  monopulse  arises 
from  the  fact  that  the  array  patterns,  particularly  in  the  radi¬ 
al  cut,  are  asymmetric.  As  a  result,  it  is  impossible  to  merely 
squint  the  two  beams  equal  amounts  from  the  nominal  array  bear¬ 
ing  and  define  the  event  location  to  be  that  angle  midway  be¬ 
tween  the  two  beams.  For  the  period  covered  by  this  report, 
asymmetric  array  patterns  are  being  handled  on  a  manual  basis; 
however,  a  computer  program  has  been  written  to  automatically 
provide  the  proper  squint  angles  for  the  computer  beamsplitting 
process . 

Finally,  there  are  two  additional  ways  in  which  the  seismic 
monopulse  system  differs  from  the  simpler  radar  one.  First, 
and  perhaps  most  important,  in  addition  to  travel-time  anomalies 
at  the  various  elements  of  1ASA,  there  are  also  amplitude  anoma¬ 
lies.  Thus,  the  signal-to-noise  ratio  at  the  various  elements 
of  the  array  may  differ  by  as  much  as  3  or  4  dB.  When  attempt¬ 
ing  to  achieve  the  maximum  performance  with  small  events,  the 
outputs  of  the  various  elements  would  have  to  be  weighted  ac¬ 
cording  to  the  received  signal-to-noise  ratios.  This  would  pro¬ 
duce  changes  in  the  basic  pattern  of  the  array  which  would  have 
to  be  accounted  for  in  the  beamsplitting  process.  Next,  it 
must  be  noted  that  all  of  the  foregoing  discussion  was  based  on 
the  fact  that  all  events  of  interest  took  place  at  the  surface 
of  the  earth.  This  is  clearly  untrue,  and  some  provision  will 
have  to  be  made  so  as  to  be  able  to  deal  with  events  occurring 
below  the  surface. 


4.3  ACCURACY  VS.  PRECISION  IN  EPI CENTRA!  DETERMINATION 

One  theoretical  point  remains  to  be  discussed  before 
examination  of  experimental  results .  This  point  has  to  do  with 
the  fact  that  the  epicenters  determined  from  a  single  1ASA  are 
subject  to  gross  errors  in  actual  location  arising  from  uncer¬ 
tainties  in  bulk  in  the  travel  times.  That  is,  while  the  pre¬ 
cision  with  which  the  event  is  located  might  become  perfect  in 
the  absence  of  noise  and  all  of  the  other  interferences  dis¬ 
cussed  before,  the  actual  epicenter  of  the  event  might  lie  in 
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some  entirely  different  place  than  that  determined  from  the 
beamsplitting  process.  As  calibration  of  the  bulk  travel  times 
of  the  earth  becomes  more  precise,  this  difficulty  may  become 
negligible . 


4*4  EXPERIMENTAL  RESULTS 

Luring  this  period,  a  computer  program  was  employed  which 
used  the  center  elements  of  the  21  LASA  clusters  to  study  three 
events.  The  first  of  these,  for  the  obvious  reason  that  it  was 
the  only  event  with  precise  position  information,  was  LONG-SHOT 
(29  Oct.  '65;  PLE  Card  No.  89;  Epicenter:  51. 4N,  179. 2E; 
Magnitude:  6.1;  Time  21  00  00.1  GMT).  The  remaining  two 
events  were  a  Kurile  Island  earthquake  (12  Lee.  *65;  PLE  Card 
No.  91;  Epicenter:  51. 5N.  178. 9W;  Lepth:  50  km;  Magnitude: 
5.0;  Time  00  48  01.7  GMT)  and  a  Semipalatinsk  surface-focus 
event(13  Eeb.  '66;  PLE  Card  No.  9;  Epicenter:  49. 8N,  78. IE; 
Magnitude:  6.3;  Time:  0  4  57.7  GMT; .  At  the  beginning  of 
this  work  it  was  not  realized  that  considerable  asymmetry 
existed  for  the  main  lobe  of  the  LASA  beam  --  particularly  in 
the  radial  direction.  For  this  reason,  the  two  beams  of  the 
monopulse  system  were  squinted  equally.  That  is,  for  example, 
in  the  azimuthal  scan  one  beam  was  squinted  6°  to  the  right, 
while  the  second  was  squinted  6°  to  the  left.  Radial  squints 
for  these  runs  were  3°.  In  addition,  the  difference  beam  out¬ 
put  was  normalized  by  dividing  by  the  unsquinted  beam  response. 
Both  of  these  procedures  were  found  to  cause  troubles  —  shifts 
in  the  null  position  of  the  monopulse  beam  and  distortion  of 
the  difference  beam  pattern.  While  it  was  not  considered  nec¬ 
essary  to  change  the  computer  program  to  take  care  of  these 
problems  during  this  period  (asymmetries  were  estimated  man¬ 
ually),  the  program  has  at  this  time  been  modified  to  auto¬ 
matically  provide  the  proper  asymmetrical  squint  by  computing 
the  array  patterns  in  the  vicinity  of  the  event  location. 

Figures  20  and  21  show  the  processed  array  outputs  while 
scanning  past  the  approximate  position  of  LONGSHOT  for  the 
azimuthal  and  radial  cases,  respectively.  As  stated  above, 
equal  azimuthal  squints  of  6°  and  radial  squint  of  3°  were 
used  for  this  computation.  The  J  and  B  surface-focus  tables 
were  used  for  these  computations.  Total  scan  for  both  cases 
includes  400  steps,  each  step  being  two  kilometers  in  range  and 
0.029°  in  bearing.  Local  travel-time  anomalies  for  this  mea¬ 
surement  were  determined  from  the  mean  value  of  the  anomalies 
measured  for  each  cluster  by  Teledyne  [11]  .  The  LONGSHOT  posi¬ 
tion  determined  from  the  zero- crossings  shown  in  these  figures 
differs  markedly  from  the  true  position  of  the  event,  approxi¬ 
mately  70  km  in  range  and  approximately  1.4°  in  bearing.  In 
addition,  the  beam  response  is  markedly  jagged  —  far  too  much 
so  for  the  signal-to-noise  ratio  of  the  received  signals. 

This  latter  phenomenon  has  since  been  traced  to  a  time-quanti¬ 
zation  noise  (time  samples  too  far  apart  for  the  differential 
type  of  calculations  being  used)  and  will  be  corrected  for  in 
the  future  by  interpolating  the  seismograms  between  samples 
as  a  part  of  the  beamsplitting  process. 
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FIGURE  20 

AZIMUTH  DIFFERENCE  PATTERN (L0NGSH0T) 


FIGURE  21 

RADIAL  DIFFERENCE  PATTERN  (10NGSH0T) 
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In  an  effort  to  discover  the  cause  of  the  location  errors 
of  this  first  calculation,  other  runs  of  the  10NGSH0T  data  were 
made.  Figures  22  through  25  give  the  outputs  of  one  of  these 
runs.  The  size  of  the  scanning  steps  for  this  rnn  are  sub¬ 
stantially  larger  than  in  the  previous  run  (0.145°  in  azimuth 
as  compared  with  0.029°,  and  0.180°  radially  as  opposed  to 
0.018°)  to  observe  gross  characteristics  of  the  array  patterns. 
Figures  26  through  29  correspond  to  tue  13  Feb.  1966  event. 

Figure  26  gives  the  unsquinted  azimuthal  scan  of  the  array, 

Figure  27,  the  difference  beam  output  for  azimuthal  scan,  Figure 
28,  the  unsquinted  radial  scan  of  the  array,  and  Figure  29,  the 
difference  beam  output  for  this  same  scan.  One  other  difference 
in  the  computation  must  be  noted.  Azimuthal  squint  for  this 
run  was  reduced  to  3°.  The  most  obvious  feature  of  this  run  is 
the  asymmetric  character  of  the  radial  unsquinted  beam.  A 
second  prominent  feature  of  this  same  pattern  is  the  existence 
of  a  substantial  flat  in  the  pattern  occurring  in  the  vicinity 
of  the  short  range  3  dB  point. 

This  asymmetry  in  pattern  coupled  with  equal  squint  angles 
will  cause  substantial  offset  in  the  zero  crossing  location  for 
the  difference  pattern  for  the  radial  case.  In  fact,  when  a 
manual  correction  for  this  asymmetry  was  applied  to  the  data 
shown  for  the  preceding  run,  the  70  km  discrepancy  between 
actual  and  measured  range  for  10NGSH0T  was  reduced  to  approxi¬ 
mately  2  km.  The  flat  in  the  pattern,  due  in  part  to  a  rela¬ 
tive  flat  in  the  travel-time  tables  for  these  ranges,  suggests 
that  extreme  care  must  be  taken  in  forming  difference  patterns. 
While  the  obvious  correction  for  such  conditions  is  to  employ 
small  squint  angles,  such  action  will  reduce  the  slope  of  the 
difference  beam  in  the  vicinity  of  the  zero  crossing,  and  hence 
lower  the  precision  of  the  measurement.  While  the  asymmetry  of 
the  unsquinted  azimuthal  beam  is  far  less  than  that  of  the  radial 
beam,  reduction  of  the  squint  angles  for  the  difference  beam 
from  six  to  three  degrees  halved  the  discrepancy  in  bearings  be¬ 
tween  actual  and  measured.  Later  refinement  of  local  travel¬ 
time  anomalies  reduced  this  error  to  less  than  0.1°  (using  ano¬ 
malies  measured  from  the  10NGSH0T  event  itself  rather  than  the 
average  anomalies  discussed  above). 

Figures  30  through  33  show  the  results  of  one  run  with  the 
small  scan  steps  for  the  earthquake.  These  runs  were  made  using 
average  anomalies  from  the  Teledyne  report.  As  can  be  seen  from 
the  figures,  the  results  of  these  computations  differ  widely 
from  the  10NGSH0T  results.  While  exceedingly  poor  azimuthal 
estimates  are  possible  from  these  data,  no  range  estimates  to 
the  events  are  possible  at  all.  It  was  suspected  that  incorrect 
travel-time  anomalies  were  responsible  for  the  poor  results 
obtained.  As  a  check  of  this  assumption,  the  Semipalatinsk 
event  was  rerun  with  the  larger  scan  steps.  Results  of  these 
computations  are  shown  in  Figures  34  through  37.  As  can  easily 
be  seen  from  these  figures,  no  real  beams  are  being  formed. 

Since  the  ability  to  form  beams  with  an  array  depends  criti¬ 
cally  on  good  timing  information,  this  lack  of  proper  beam  for¬ 
mation  suggests  strongly  that  the  travel-time  anomalies  were 
in  fact  incorrect. 
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The  results  outlined  above  suggest  the  directions  to  be 
taken  in  the  next  period  for  the  problem  of  beam  splitting. 
First,  more  accurate  local  travel-time  anomalies  must  be  sought 
This  will  not  only  include  anomalies  for  the  events  studied  so 
far,  but  also  average  anomalies  for  several  regions  under 
study.  Second,  time  interpolation  of  seismograms  will  be  in¬ 
cluded  in  the  computer  program  to  minimize  the  time-quantiza¬ 
tion  noise  described  above.  Third,  automatic  asymmetric  squint 
ing  will  be  included  in  the  program.  Fourth,  normalization  of 
the  difference  beam  by  the  unsquinted  pattern  of  the  array  will 
be  eliminated.  Finally,  polynomial  fitting  of  the  difference 
beam  to  improve  estimates  of  the  zero  crossings  will  also  be 
included  in  the  program. 

In  addition  to  these  corrections  to  the  computational 
procedure,  the  number  of  events  studied  will  be  increased  with 
emphasis  placed  on  small  magnitude  events. 
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FIGURE  22 
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UNSQUINTED  RADIAL  PATTERN  (LONGSHOT) 


RADIAL  DIPFERENCE  PATTERN  (LONGSHOT) 
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FIGURE  26 

UNSQUINTED  AZIMUTH  PATTERN  (SEMIPALATINSK  13  FEB.  ’66) 


AZIMUTH  DIFFERENCE  PATTERN  (SEMIPALATINSK  13  FEB.  '66) 
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UNSQUINTED  RADIAL  PATTERN  ( S EMI PALATINS K  13  EEB.  *66) 


FIGURE  29 

RADIAL  DIFFERENCE  PATTERN  (SEMIPALATINSK  13  FEB.  *66) 
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UNSQUINTED  AZIMUTH  PATTERN  (E.Q .  KURILE  IS.  12  DEC.  *65) 
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FIGURE  32 

UNSQUINTED  RADIAL  PATTERN  (E.Q.  KURILE  IS.  12  DEC,  '65) 


FIGURE  33 

RADIAL  DIFFERENCE  PATTERN  (E.Q.  KURILE  IS.  12  DEC.  ’65) 


FIGURE  34 

UNSQUINTED  AZIMUTH  PATTERN  (SEMIPA1ATINSK,  FEB  13  '66) 
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UNSQUINTED  RADIAL  PATTERN  ( S EMI PAL ATI NS K,  PEB  13  '66) 


j _ 1  FIGURE  37 _ , 

RADIAL  DIFFERENCE  PATTERN  (SMIPALATINSK,  FEB  13  *66) 
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